// John Lake
// CSCE 441

#include <stdlib.h>
#include <math.h>
#include "RayTracer.h"

Point EyeLocation = Point(ImageW/2, ImageH/2, -ImageD);

/*************************************************************************
 Applies the lighting equation given a normal vector.  In this assignment
   all the coefficients and light sources are given and thus explicitly 
   defined here, though this could be generalized to accept coefficients
   and light sources

Does NOT do the ambient lighting

 p : Point at which the surface is intersected
 n : normal to the surface at p

 LightingPackage contains the following:
 lights : vector of the light sources
 ambientLight : the ambient light present in the scene
 objects : vector of objects present in the scene (used for shadow calculations)
				(must be pointers since is virtual class)
 coeffs : containing ka, kd, ks (in that order) for the object
*************************************************************************/
Color lightingEquation(Point p, Normal n, Color kd, Color ks, LightingPackage lp)
{
	std::vector<Light> lights = lp.lights;	
	std::vector<Object*> objects = lp.objects;
	Color diffuse, specular, final;
	
	// clear the colors to begin
	diffuse = Color();
	specular = Color();
	final = Color();


	std::vector<Light>::iterator lIter;
	std::vector<Object*>::iterator oIter;

	Light currentLight;


	for(lIter = lights.begin(); lIter != lights.end(); lIter++)
	{
		currentLight = (*lIter);

		// find L : direction to light
		Normal lightDirection = currentLight.location() - p;
		lightDirection.normalize();

		Ray toLight = Ray(p, lightDirection);
		
		// find E : direction to viewer
		Normal eye = EyeLocation - p;
		eye.normalize();

		bool flag = FALSE;


		// check shadows against location of light (else infinite planes will always intersect
		float lightParam = (*lIter).intersection(toLight);

		// check for shadows
		for(oIter = objects.begin(); oIter != objects.end(); oIter++)
		{
			if(flag == TRUE)
				break;

			Ray ray = Ray(p, (*lIter).location() - p);	// draw a Ray from the location to the light, see if there's an intersection for each object
			Ray invDefRay = Ray((*oIter)->deformationInverse*ray.p(), (*oIter)->deformationInverse*ray.v());

			float check = (*oIter)->intersection(invDefRay);


			// strange number below because of problems with non-exact numbers
			if(check >= TOLERANCE)	// there was an intersection => in shadow for this light source
			{
				if(!(lightParam <= check))	// if is a shorter distance to the light than the obstruction
					flag = TRUE;
			}
		}

		if(flag == TRUE)	// if in shadow for this light source, no need to do other equations
			continue;


		// Diffuse
		float LdotN = lightDirection*n;
		diffuse.r = currentLight.color().r*kd.r*LdotN;
		diffuse.g = currentLight.color().g*kd.g*LdotN;
		diffuse.b = currentLight.color().b*kd.b*LdotN;
		diffuse.clamp();


		// Specular
			// R = 2*(L.N)N - L
		Normal sRefl;	// reflected specular light
		sRefl.x = 2*LdotN*n.x - lightDirection.x;
		sRefl.y = 2*LdotN*n.y - lightDirection.y;
		sRefl.z = 2*LdotN*n.z - lightDirection.z;
		sRefl.normalize();

		float RdotE = sRefl * eye;
		float alpha = pow(RdotE,SPEC_EXP);	// raise R.E to the specular coefficient's power
			// specular exponent is always odd, to preserve sign

		specular.r = currentLight.color().r*ks.r*alpha;
		specular.g = currentLight.color().g*ks.g*alpha;
		specular.b = currentLight.color().b*ks.b*alpha;
		specular.clamp();
	
		final.r += diffuse.r + specular.r;
		final.g += diffuse.g + specular.g;
		final.b += diffuse.b + specular.b;
	}


	final.clamp();

	return final;


}

// Note - this converts to spherical as if the z direction was in what is really the y direction
//   thus, it looks more natural to us
SpherePoint toSpherical(Point originalPoint, Point center)
{
	// radius is an extraneous argument
	SpherePoint rvalue;

	Point p = originalPoint - center;	// translate so that the center is (0,0,0)
	
	rvalue.center = center;
	rvalue.radius = sqrt(p.x*p.x + p.y*p.y + p.z*p.z);
	if(p.x == p.z)
		rvalue.theta = 0;
	else
		rvalue.theta = atan(p.x/p.z);
	rvalue.phi = acos(p.y/rvalue.radius);
	return rvalue;

}




/***********************************  Color  ******************************************/
// Clamps RGB values to [0,1]
void 
Color::clamp()
{
	if(r < 0)
		r = 0;
	if(g < 0)
		g = 0;
	if(b < 0)
		b = 0;
	
	if(r > 1)
		r = 1;
	if(g > 1)
		g = 1;
	if(b > 1)
		b = 1;
}


Color 
Color::operator+(Color& rhs)
{
	Color rvalue;
	rvalue.r = this->r + rhs.r;
	rvalue.g = this->g + rhs.g;
	rvalue.b = this->b + rhs.b;
	return rvalue;
}


Color 
Color::operator*(float& rhs)
{
	Color rvalue;
	rvalue.r = this->r * rhs;
	rvalue.g = this->g * rhs;
	rvalue.b = this->b * rhs;
	return rvalue;
}



/***********************************  Point  ******************************************/

// Normalizes a vector; not really used in a purely "Point" context
void
Point::normalize()
{
	float magnitude;
	
	magnitude = sqrt(x*x + y*y + z*z);
	
	x /= magnitude;
	y /= magnitude;
	z /= magnitude;
}


// get the magnitude of a vector
float
Point::magnitude()
{
	return sqrt(x*x + y*y + z*z);
}


/**********************************************************
  Does a component-wise subtraction, useful for drawing a
   vector between points
/**********************************************************/
Point 
Point::operator-(Point& rhs)
{
	Normal rvalue;
	rvalue.x = this->x - rhs.x;
	rvalue.y = this->y - rhs.y;
	rvalue.z = this->z - rhs.z;
	return rvalue;
}


/**********************************************************
 Does a component-wise addition 
/**********************************************************/
Point 
Point::operator+(Point& rhs)
{
	Normal rvalue;
	rvalue.x = this->x + rhs.x;
	rvalue.y = this->y + rhs.y;
	rvalue.z = this->z + rhs.z;
	return rvalue;
}


/**********************************************************
 Does dot multiplication
 /**********************************************************/
float 
Point::operator*(Point& rhs)
{
	float rvalue;
	rvalue = this->x * rhs.x + this->y * rhs.y + this->z * rhs.z;
	return rvalue;
}


/**********************************************************
 Does scalar multiplication
/**********************************************************/
Point 
Point::operator*(float& rhs)
{
	Normal rvalue;
	rvalue.x = this->x * rhs;
	rvalue.y = this->y * rhs;
	rvalue.z = this->z * rhs;
	return rvalue;
}


/**********************************************************
 Does scalar divide
/**********************************************************/
Point 
Point::operator/(float& rhs)
{
	Point rvalue;
	rvalue.x = this->x / rhs;
	rvalue.y = this->y / rhs;
	rvalue.z = this->z / rhs;
	return rvalue;
}


/***********************************  Ray  ******************************************/
Point
Ray::p()
{
	return point;
}


Normal
Ray::v()
{
	return vector;
}


Point
Ray::at(float parameter)
{
	Point rvalue;
	rvalue = point + vector*parameter;
	return rvalue;
}


float
Ray::magnitude()
{
	return vector.magnitude();
}


/***********************************  Light  ******************************************/
Point 
Light::location()
{  return l;  }


Color 
Light::color()
{  return c;  }


// This is a really roundabout way to do this, but it works
float 
Light::intersection(Ray ray)
{
	float t;
	
	float alpha = ray.v() * (l - ray.p());
	t = alpha / (ray.v() * ray.v());

	return t;
}


/***********************************  Matrix  ******************************************/

Matrix 
Matrix::transpose()
{
	Matrix rvalue;

	rvalue.data[0][0] = data[0][0];
	rvalue.data[0][1] = data[1][0];
	rvalue.data[0][2] = data[2][0];

	rvalue.data[1][0] = data[0][1];
	rvalue.data[1][1] = data[1][1];
	rvalue.data[1][2] = data[2][1];
	
	rvalue.data[2][0] = data[0][2];
	rvalue.data[2][1] = data[1][2];
	rvalue.data[2][2] = data[2][2];
	
	return rvalue;
}


Point 
Matrix::operator*(Point& rhs)
{
	Point rvalue;

	rvalue.x = data[0][0]*rhs.x + data[0][1]*rhs.y + data[0][2]*rhs.z;
	rvalue.y = data[1][0]*rhs.x + data[1][1]*rhs.y + data[1][2]*rhs.z;
	rvalue.z = data[2][0]*rhs.x + data[2][1]*rhs.y + data[2][2]*rhs.z;

	return rvalue;
}


/***********************************  Object  ******************************************/
Object::Object()
{
	// set the transform to I
	deformationInverse.data[0][0] = 1;
	deformationInverse.data[1][1] = 1;
	deformationInverse.data[2][2] = 1;
	
	// set the coefficients to 0 to begin with - user should change
	ka = kd = ks = Color(0,0,0);

	refl_coeff = 0;
}


void 
Object::setDeformationInverse(Matrix _deformationInverse)
{
	deformationInverse = _deformationInverse;
}


void 
Object::setKa(Color _ka)
{  ka = _ka;  }

void 
Object::setKd(Color _kd)
{  kd = _kd;  }

void 
Object::setKs(Color _ks)
{  ks = _ks;  }

void
Object::setReflCoeff(float _gamma)
{  refl_coeff = _gamma;  }


Color 
Object::getKa(Point p)
{  return ka;  }

Color
Object::getKd(Point p)
{  return kd;  }

Color 
Object::getKs(Point p)
{  return ks;  }

Color
Object::lightingFunction(Ray ray, float t, float depth, Color ambientLight, LightingPackage lp)
{
	std::vector<Object*> objects = lp.objects;
	Color final = Color(0,0,0);

	// 0. Check depth
	if(depth > MAX_DEPTH)
		return Color(0,0,0);

	// 1. Retrieve normal
	Ray invDefRay = Ray(deformationInverse*ray.p(), deformationInverse*ray.v());
	Point p = invDefRay.at(t);	// pretty sure I need to do it like this to put it in the non-deformed space
	Normal n = getNormal(p);	// inverse transpose stuff done inside getNormal
	n.normalize();

	// 2. Retrieve ka, ks, kd
	Color ka = getKa(p);
	Color kd = getKd(p);
	Color ks = getKs(p);

	// 3. Apply lighting equation
	// Lighting Equation: I = ka*A + C(kd(L.N)+ks(R.E)^n)

	//  a. do ambient lighting
		// Ambient light - independent of light source
	Color ambient;
	ambient.r = ka.r*ambientLight.r;
	ambient.g = ka.g*ambientLight.g;
	ambient.b = ka.b*ambientLight.b;
	ambient.clamp();

	Color diffuseAndSpecular = lightingEquation(p, n, kd, ks, lp);
	diffuseAndSpecular.clamp();

	final = ambient + diffuseAndSpecular;	

	// 4. recurse/reflect
	//  a. find reflected ray
	float beta = (ray.v() * n) * 2;
	Normal gamma = n * beta;
	Normal reflected = ray.v() - gamma;	// negative of powerpoint slides
	
	Ray reflectedRay = Ray(p, reflected);

	//  b. intersect?
	std::vector<Object*>::iterator oIter;
	int indexofIntersected = -1;
	int parameterofIntersected = -1;

	for(oIter = objects.begin(); oIter != objects.end(); oIter++)
	{
		// no self (i.e. internal) reflections
		if((*oIter)->index == index)
			continue;

		Ray invDefReflRay = Ray((*oIter)->deformationInverse*reflectedRay.p(), (*oIter)->deformationInverse*reflectedRay.v());
		float parameter = (*oIter)->intersection(invDefReflRay);
		if(parameter >= TOLERANCE)	// intersection!
		{
			if((parameter < parameterofIntersected) || (parameterofIntersected == -1))	// closer (or new) object
			{
				parameterofIntersected = parameter;
				indexofIntersected = (*oIter)->index;
			}
		}
	}

	if(indexofIntersected != -1)
	{	// reflection intersections!
		for(oIter = objects.begin(); oIter != objects.end(); oIter++)
			if((*oIter)->index == indexofIntersected)
				break;

		//  c. recurse (increasing depth)
		final = final + (*oIter)->lightingFunction(reflectedRay, parameterofIntersected, depth+1, ambientLight, lp) * refl_coeff;	// scalar multiply
		final.clamp();
	}

	// 5. return final?
	return final;
}

/***********************************  Plane  ******************************************/
float 
Plane::intersection(Ray ray)
{
	float t;
	float nDotV = normal*ray.v();
	
	if (nDotV == 0) {
		// Ray is parallel to plane - no intersection
		return -1;
	}
	
	// there is only one pointOfIntersection for planes
	t = normal * (pointOnPlane - ray.p())/nDotV;	

	return t;
}

Normal 
Plane::getNormal(Point p)
{
	Normal rvalue = deformationInverse.transpose() * normal;
	return rvalue;  
}


/***********************************  Cylinder  ******************************************/
float 
Cylinder::intersection(Ray ray)
{
	// First, project everything onto a plane
	float pdota = ray.p() * axis;	// C++ angry about operands again	
	float vdota = ray.v() * axis;
	Point pplane = ray.p() - axis * pdota;		//(ray.p() * axis);
	Normal vplane = ray.v() - axis * vdota;		//(ray.v() * axis);
	vplane.normalize();
	
	Ray rplane = Ray(pplane, vplane);	// ray in the plane orthognal to the axis

	float cdota = center * axis;
	Point cplane = center - axis * cdota;	// center point on the projected plane


	// greek letters because C++ was getting angry at my operators
	Normal beta = cplane - rplane.p();	// component-wise subtraction	// original
	float gamma = beta * rplane.v();	// dot product
	float delta = rplane.v() * rplane.v();	// dot product
	Normal epsilon = rplane.v() * gamma;	// scalar product;
	Normal W = rplane.p() - cplane + (rplane.v() * gamma) / delta;	// component-wise addition, subtraction; scalar multiplication, division
//	Normal W = ray.p() - center + (ray.v() * ((center - ray.p()) * ray.v())) / (ray.v() * ray.v());

	float t1, t2, alpha;

	if (W.magnitude() > radius) {
		// no intersections
		return -1;
	}
	else if (W.magnitude() == radius) {
		// 1 intersection - tangential
		float t1 = (center - ray.p())*ray.v()/(ray.v()*ray.v());

		return t1;
	}
	else if (W.magnitude() < radius) {
		// 2 intersections - find closer intersection
		alpha = sqrt(pow(radius,2) - pow(W.magnitude(),2));
		
		t1 = alpha/rplane.magnitude() + (cplane - rplane.p())*rplane.v()/(rplane.v()*rplane.v());
		
		t2 = - alpha/rplane.magnitude() + (cplane - rplane.p())*rplane.v()/(rplane.v()*rplane.v());

		if(t1 == t2)
		{
			return -1;	// should never happen
		}

		if(t1 < t2)
			return t1;
		else 
			return t2;

	}
	else
	{
		return -1;	// should never happen
	}


}

Normal 
Cylinder::getNormal(Point p)
{  
	Normal rvalue;
	Normal beta = p - center;
	float gamma = beta * axis;
	Normal delta = axis * gamma;
	Normal temp = beta - delta;	// = (p-center) - axis * (p-center * axis);
	rvalue = deformationInverse.transpose() * temp;
	return rvalue;
}


/***********************************  Sphere  ******************************************/
float 
Sphere::intersection(Ray ray)
{
	// greek letters because C++ was getting angry at my operators
	Normal beta = center - ray.p();	// component-wise subtraction
	float gamma = beta * ray.v();	// dot product
	float delta = ray.v() * ray.v();	// dot product
	Normal epsilon = ray.v() * gamma;	// scalar product;
	Normal W = ray.p() - center + (ray.v() * gamma) / delta;	// component-wise addition, subtraction; scalar multiplication, division
//	Normal W = ray.p() - center + (ray.v() * ((center - ray.p()) * ray.v())) / (ray.v() * ray.v());
	
	if (W.magnitude() > radius) {
		// no intersections
		return -1;
	}
	else if (W.magnitude() == radius) {
		// 1 intersection - tangential
		float t = (center - ray.p())*ray.v()/(ray.v()*ray.v());

		return t;
	}
	else if (W.magnitude() < radius) {
		// 2 intersections - find closer intersection
		float alpha = sqrt(pow(radius,2) - pow(W.magnitude(),2));
		
		float t1 = alpha/ray.magnitude() + (center - ray.p())*ray.v()/(ray.v()*ray.v());
		
		float t2 = -alpha/ray.magnitude() + (center - ray.p())*ray.v()/(ray.v()*ray.v());
		
		if(t1 == t2)
		{
			return -1;	// should never happen
		}

		if(t1 < t2)
			return t1;
		else 
			return t2;

	}
	else
	{
		return -1;	// should never happen
	}
	
}

Normal 
Sphere::getNormal(Point p)
{
	Point temp = p - center;
	Point rvalue = deformationInverse.transpose() * temp;
	return rvalue;
}


/***********************************  Earth  ******************************************/
bool
Earth::isOnContinent(SpherePoint sp)
{
	// The continents are effectively a pattern of squares

	// I was having problems with the logic, else I would have returned true
	for(int i = -16; i < 16; i += 2)
		if(((sp.theta > i*PI/16) && (sp.theta < (i+1)*PI/16)) || ((sp.phi > i*PI/16) && (sp.phi < (i+1)*PI/16)))
			return false;

	return true;
}


bool
Earth::isOnIceCap(SpherePoint sp)
{
	if((sp.phi < PI/5) || (sp.phi > 4*PI/5))
		return true;

	return false;
}


Color 
Earth::getKa(Point p)
{
	return Color(0.1,0.1,0.1);
}


Color 
Earth::getKd(Point p)
{
	SpherePoint loc = toSpherical(p,center);
	if(isOnIceCap(loc))
		return Color(0.8,0.8,0.8);

	srand(p.x+p.y+p.z);	// gonna want the colors to be consistent for reflections and stuff
	float a = rand();
	float b = rand();

	bool Clouds;

	if(a/b + 0.9 < 1)	// sea
		Clouds = false;
	else		// clouds
		Clouds = true;

	if(isOnContinent(loc))
		if(Clouds)
			return Color(0.6,0.7,0.6);
		else
			return Color(0,0.7,0.1);
	else
		if(Clouds)
			return Color(0.7,0.7,0.7);
		else
			return Color(0,0.1,0.7);

}


Color 
Earth::getKs(Point p)
{
	SpherePoint loc = toSpherical(p,center);

	if(isOnContinent(loc))
		return Color(0,0,0);
	else return Color(0.6,0.6,0.6);
}


/***********************************  Jupiter  ******************************************/


bool 
Jupiter::isRedSpot(Point p)
{
	Point spotCenter = center + Point(radius/3.5,-radius/4,radius/4);
	float spotRadius = radius/10;

	Point test = p - spotCenter;	// gives the distance from the center of the point

	if(test.x*test.x/4 + test.y*test.y < spotRadius*spotRadius)	// inside an ellipse
		return true;

	return false;
	
}


bool 
Jupiter::isYellow(SpherePoint sp)
{
	for(int i = 0; i < 40; i += 4)
		if((sp.phi > i*PI/40) && (sp.phi < (i+1)*PI/40))
			return true;

	return false;

}


bool 
Jupiter::isOrange(SpherePoint sp)
{
	for(int i = 1; i < 40; i += 4)
		if((sp.phi > i*PI/40) && (sp.phi < (i+1)*PI/40))
			return true;

	return false;
}


bool 
Jupiter::isWhite(SpherePoint sp)
{
	for(int i = 2; i < 40; i += 4)
		if((sp.phi > i*PI/40) && (sp.phi < (i+1)*PI/40))
			return true;

	return false;
}


bool 
Jupiter::isBrown(SpherePoint sp)
{
	for(int i = 3; i < 40; i += 4)
		if((sp.phi > i*PI/40) && (sp.phi < (i+1)*PI/40))
			return true;

	return false;
}

Color 
Jupiter::getKa(Point p)
{
	return Color(0.1,0.1,0.1);
}

Color 
Jupiter::getKd(Point p)
{
	SpherePoint loc = toSpherical(p,center);
	if(isRedSpot(p))
		return Color(0.8,0.1,0);

	if(isYellow(loc))
		return Color(0.8,0.7,0);

	if(isOrange(loc))
		return Color(0.8,0.4,0);

	if(isWhite(loc))
		return Color(0.7,0.6,0.6);

	if(isBrown(loc))
		return Color(0.4,0.2,0);

	
	return Color(0,0,0);
}


Color 
Jupiter::getKs(Point p)
{
	return Color(0,0,0);
}


/***********************************  Neptune  ******************************************/

bool 
Neptune::isDarkSpot(Point p)
{
	Point spotCenter = center + Point(radius/3.5,-radius/4,radius/4);
	float spotRadius = radius/10;

	Point test = p - spotCenter;	// gives the distance from the center of the point

	if(test.x*test.x/4 + test.y*test.y < spotRadius*spotRadius)	// inside an ellipse
		return true;

	return false;
	
}


Color 
Neptune::getKa(Point p)
{
	return Color(0.05,0.05,0.05);
}

Color 
Neptune::getKd(Point p)
{
	if(isDarkSpot(p))
		return Color(0,0,0.7);
	
	return Color(0,0,0.9);
}


Color 
Neptune::getKs(Point p)
{
	return Color(0,0,0);
}


/***********************************  Ring  ******************************************/
Ring::Ring(Point _center, float _radius1, float _radius2, Normal _normal, int _index)
{
	inner = Sphere(_center, _radius1, 1);
	outer = Sphere(_center, _radius2, 2);
	plane = Plane(_normal, _center, 3);

	index = _index;
}


Normal 
Ring::getNormal(Point p)
{
	return plane.getNormal(p);
}


float 
Ring::intersection(Ray ray)
{
	float tPlane = plane.intersection(ray);	// gives the parameter and consequent location of the intersection

	Point p = ray.at(tPlane);	// Point of intersection

	SpherePoint sp = toSpherical(p,inner.center);	// whether the inner or outer's center is used is irrelevant

	if(sp.radius > outer.radius)
	{	// no intersection - too far out
		return -1;
	}
	else	// sp.radius < outer.radius
	{
		if(sp.radius < inner.radius)
		{	// no intersection - too close to the center
			return -1;
		}
		else
		{	// is between the first and second radii - ring present!
			return tPlane;
		}
	}
}
